##################################################
#
# PLOTS OF MEASURES OF WOMEN'S STATUS/RIGHTS
# MEASURED AS CONTINUOUS VARIABLES FOR THE 
# OIL, ISLAM, and WOMEN PROJECT
#
# Paul Musgrave and Yu-Ming Liou
# musgrave@umass.edu and yl254@georgetown.edu
#
# Created 27 February 2016
#
##################################################


##################################################
# R Housekeeping
##################################################

# Remove and empty workspace
rm(list=ls())

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Packages
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

library(foreign)
library(sandwich)
library(apsrtable)
library(tonymisc)  # includes summaryR()
library(MASS) # for polr() which runs our ologit

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Functions
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 



# Change this to  your working directory
mywd <- "/Users/paulmusgrave/Dropbox/0001 Academic Projects/Ongoing/0127 Oil Islam Women/ISQ Accepted Submission/Replication"

setwd(paste(mywd,"/Code",sep=""))

source("statifyFN.R")
source("multipleOLSFN.R")
source("multipleOlogitFN.R")
source("genModelsFN.R")
source("plotlineFN.R")
source("plotprepFN.R")
source("modelDefinitions.R")

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Load Data
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 

setwd(paste(mywd,"/Data",sep=""))

data <- read.dta("NewDVModelsXSforISQFINAL2015.dta")


# NOTE THAT FACTORS HAVE BEEN REORDERED SO THAT
# WORSE OUTCOMES FOR WOMEN ARE CONSISTENTLY ** LOWER **

data$ciri_wecon_02 <- factor(data$ciri_wecon_02,
                             levels=c(0,1,2,3),
                             ordered=TRUE)
data$ciri_wopol_02 <- factor(data$ciri_wopol_02,
                             levels=c(0,1,2,3),
                             ordered=TRUE)
data$ciri_wosoc_02 <- factor(data$ciri_wosoc_02,
                             levels=c(0,1,2,3),
                             ordered=TRUE)

data$ciri_wecon_04 <- factor(data$ciri_wecon_04,
                             levels=c(0,1,2,3),
                             ordered=TRUE)
data$ciri_wopol_04 <- factor(data$ciri_wopol_04,
                             levels=c(0,1,2,3),
                             ordered=TRUE)
data$ciri_wosoc_04 <- factor(data$ciri_wosoc_04,
                             levels=c(0,1,2,3),
                             ordered=TRUE)

data$McDermott_07 <- factor(data$McDermott_07,
#                           levels=c(0,1,2,3,4),
                            levels=c(4:0),
                            ordered=TRUE)
data$Caprioli_07 <- factor(data$Caprioli_07,
#                           levels=c(0,1,2,3,4),
                           levels=c(4:0),
                           ordered=TRUE)

data$McDermott_11 <- factor(data$McDermott_11,
 #                           levels=c(0,1,2,3,4),
                            levels=c(4:0),
                            ordered=TRUE)
data$Brinton_10 <- factor(data$Brinton_10,
#                          levels=c(1:4),
                          levels=c(4:1),
                          ordered=TRUE)

data$Clan.fac<- factor(data$Clan,
#                          levels=c(1:12),
                            levels=c(12:1),
                          ordered=TRUE)

data$GII_UN_11 <- -(data$GII_UN_11) # to reorder UN so everything is negative

data$oil_gas_pc_11_log <- log(data$oil_gas_pc_11+1)
data$oil_gas_pc_07_log <- log(data$oil_gas_pc_07+1)
data$oil_gas_pc_04_log <- log(data$oil_gas_pc_04+1)
data$oil_gas_pc_02_log <- log(data$oil_gas_pc_02+1)


data$oil_gas_pc_11 <- data$oil_gas_pc_11 / 1000
data$oil_gas_pc_07 <- data$oil_gas_pc_07 / 1000
data$oil_gas_pc_04 <- data$oil_gas_pc_04 / 1000
data$oil_gas_pc_02 <- data$oil_gas_pc_02 / 1000

#################################################
#
# New DVs
# 
#################################################

# Brinton
out.Brinton.10<- multipleOlogit(formulas=genModels(c("Brinton_10"),
                                                   c(p140.11,p142.11)),
				          	        modelnames=c(LETTERS[1:7]),
                            y=data)

out.Brinton.10.log<- multipleOlogit(formulas=genModels(c("Brinton_10"),
                                                   c(p140.11.log,p142.11.log)),
                                modelnames=c(LETTERS[1:7]),
                                y=data)


# McDermott 11
out.McDermott.11<- multipleOlogit(formulas=genModels(c("McDermott_11"),
                                                     c(p140.11,p142.11)),
                                modelnames=c(LETTERS[1:7]),
                                y=data)


out.McDermott.11.log<- multipleOlogit(formulas=genModels(c("McDermott_11"),
                                                       c(p140.11.log,
                                                         p142.11.log)),
                                      modelnames=c(LETTERS[1:7]),
                                      y=data)


# CLAN done two ways
out.Clan.11<- multipleOlogit(formulas=genModels(c("Clan.fac"),
                                                c(p140.11,
                                                  p142.11)),
                             modelnames=c(LETTERS[1:7]),
                             y=data)

out.Clan.11.log<- multipleOlogit(formulas=genModels(c("Clan.fac"),
                                                c(p140.11.log,
                                                  p142.11.log)),
                                 modelnames=c(LETTERS[1:7]),
                                 y=data)


out.Clan.11.ols<- multipleOLS(formulas=genModels(c("Clan"),
                                                c(p140.11,
                                                  p142.11)),
                              modelnames=c(LETTERS[1:7]),
                              y=data)

out.Clan.11.ols.log<- multipleOLS(formulas=genModels(c("Clan"),
                                                    c(p140.11.log,
                                                      p142.11.log)),
                                  modelnames=c(LETTERS[1:7]),
                                  y=data)


# UN GII
out.GII.ols<- multipleOLS(formulas=genModels(c("GII_UN_11"),
                                                 c(p140.11,
                                                   p142.11)),
                          modelnames=c(LETTERS[1:7]),
                          y=data)

out.GII.11.ols.log<- multipleOLS(formulas=genModels(c("GII_UN_11"),
                                                     c(p140.11.log,p142.11.log)),
                                 modelnames=c(LETTERS[1:7]),
                                 y=data)

# McDermott 2007
out.McDermott_07<- multipleOlogit(formulas=genModels(c("McDermott_07"),
                                                     c(p140.07,p142.07)),
                                  modelnames=c(LETTERS[1:7]),
                                  y=data)


out.McDermott_07.log<- multipleOlogit(formulas=genModels(c("McDermott_07"),
                                                         c(p140.07.log,
                                                           p142.07.log)),
                                      modelnames=c(LETTERS[1:7]),
                                      y=data)

# Caprioli_07

out.Caprioli_07<- multipleOlogit(formulas=genModels(c("Caprioli_07"),
                                                    c(p140.07,p142.07)),
                                 modelnames=c(LETTERS[1:7]),
                                 y=data)


out.Caprioli_07.log<- multipleOlogit(formulas=genModels(c("Caprioli_07"),
                                                        c(p140.07.log,
                                                          p142.07.log)),
                                     modelnames=c(LETTERS[1:7]),
                                     y=data)


# CIRI WOSOC 04

out.ciri.wosoc.04<- multipleOlogit(formulas=genModels(c("ciri_wosoc_04"),
                                                      c(p140.04,p142.04)),
                                   modelnames=c(LETTERS[1:7]),
                                   y=data)


out.ciri.wosoc.04.log<- multipleOlogit(formulas=genModels(c("ciri_wosoc_04"),
                                                          c(p140.04.log,
                                                            p142.04.log)),
                                       modelnames=c(LETTERS[1:7]),
                                       y=data)

# CIRI WOSOC 2002
out.ciri.wosoc.02<- multipleOlogit(formulas=genModels(c("ciri_wosoc_02"),
                                                      c(p140.02,p142.02)),
                                   modelnames=c(LETTERS[1:7]),
                                   y=data)


out.ciri.wosoc.02.log<- multipleOlogit(formulas=genModels(c("ciri_wosoc_02"),
                                                          c(p140.02.log,
                                                            p142.02.log)),
                                       modelnames=c(LETTERS[1:7]),
                                       y=data)

# CIRI wopol 04

out.ciri.wopol.04<- multipleOlogit(formulas=genModels(c("ciri_wopol_04"),
                                                      c(p140.04,p142.04)),
                                   modelnames=c(LETTERS[1:7]),
                                   y=data)


out.ciri.wopol.04.log<- multipleOlogit(formulas=genModels(c("ciri_wopol_04"),
                                                          c(p140.04.log,
                                                            p142.04.log)),
                                       modelnames=c(LETTERS[1:7]),
                                       y=data)

# CIRI wopol 2002
out.ciri.wopol.02<- multipleOlogit(formulas=genModels(c("ciri_wopol_02"),
                                                      c(p140.02,p142.02)),
                                   modelnames=c(LETTERS[1:7]),
                                   y=data)


out.ciri.wopol.02.log<- multipleOlogit(formulas=genModels(c("ciri_wopol_02"),
                                                          c(p140.02.log,
                                                            p142.02.log)),
                                       modelnames=c(LETTERS[1:7]),
                                       y=data)


# CIRI wecon 04

out.ciri.wecon.04<- multipleOlogit(formulas=genModels(c("ciri_wecon_04"),
                                                      c(p140.04,p142.04)),
                                   modelnames=c(LETTERS[1:7]),
                                   y=data)


out.ciri.wecon.04.log<- multipleOlogit(formulas=genModels(c("ciri_wecon_04"),
                                                          c(p140.04.log,
                                                            p142.04.log)),
                                       modelnames=c(LETTERS[1:7]),
                                       y=data)

# CIRI wecon 2002
out.ciri.wecon.02<- multipleOlogit(formulas=genModels(c("ciri_wecon_02"),
                                                      c(p140.02,p142.02)),
                                   modelnames=c(LETTERS[1:7]),
                                   y=data)


out.ciri.wecon.02.log<- multipleOlogit(formulas=genModels(c("ciri_wecon_02"),
                                                          c(p140.02.log,
                                                            p142.02.log)),
                                       modelnames=c(LETTERS[1:7]),
                                       y=data)



## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
# Plotting Models
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #
vars.to.plot <- c("oil_gas_pc_02",
                  "oil_gas_pc_04",
                  "oil_gas_pc_07",
                  "oil_gas_pc_11")
				

plot.out.ciri.wecon.02 <- do.call(rbind.data.frame,
                             lapply(out.ciri.wecon.02,plotprep,vars.to.plot))
plot.out.ciri.wecon.04 <- do.call(rbind.data.frame,
                                 lapply(out.ciri.wecon.04,
                                        plotprep,
                                        vars.to.plot))
plot.out.ciri.wopol.04 <- do.call(rbind.data.frame,
                                 lapply(out.ciri.wopol.04,
                                        plotprep,
                                        vars.to.plot))
plot.out.ciri.wopol.02 <- do.call(rbind.data.frame,
                                 lapply(out.ciri.wopol.02,
                                        plotprep,
                                        vars.to.plot))
plot.out.ciri.wosoc.04 <- do.call(rbind.data.frame,
                                  lapply(out.ciri.wosoc.04,
                                         plotprep,
                                         vars.to.plot))
plot.out.ciri.wosoc.02 <- do.call(rbind.data.frame,
                                  lapply(out.ciri.wosoc.02,
                                         plotprep,
                                         vars.to.plot))
plot.out.Caprioli_07 <- do.call(rbind.data.frame,
                                lapply(out.Caprioli_07,
                                       plotprep,
                                       vars.to.plot))
plot.out.McDermott_07 <- do.call(rbind.data.frame,
                                lapply(out.McDermott_07,
                                       plotprep,
                                       vars.to.plot))
plot.out.GII.ols <- do.call(rbind.data.frame,
                            lapply(out.GII.ols,
                                   plotprep,
                                   vars.to.plot))
plot.out.Clan.11.ols <- do.call(rbind.data.frame,
                                lapply(out.Clan.11.ols,
                                       plotprep,
                                       vars.to.plot))
plot.out.Clan.11 <- do.call(rbind.data.frame,
                            lapply(out.Clan.11,
                                   plotprep,
                                   vars.to.plot))
plot.out.McDermott.11 <- do.call(rbind.data.frame,
                                 lapply(out.McDermott.11,
                                        plotprep,
                                        vars.to.plot))
plot.out.Brinton.10 <- do.call(rbind.data.frame,
                               lapply(out.Brinton.10,
                                      plotprep,
                                      vars.to.plot))

moreprep <- function(x){
  x$Index <- 1:nrow(x)
  x$Significance <- "white"
  x$Significance[abs(x$"T-Score")>=1.96] <- "black"
  x$Significance[abs(x$"T-Score")>=1.65 & abs(x$"T-Score")<1.96] <- "gray"
  return(x)
}


moreprepOLS <- function(x){
  x$Index <- 1:nrow(x)
  x$Significance <- "white"
  x$"T-Score" <- (x$Estimate)/(x$SE)
  x$Significance[abs( x$"T-Score")>=1.96] <- "black"
  x$Significance[abs( x$"T-Score")>=1.65 & abs( x$"T-Score")<1.96] <- "gray"
  return(x)
}

plot.out.ciri.wecon.02 <- moreprep(plot.out.ciri.wecon.02)
plot.out.ciri.wecon.04 <- moreprep(plot.out.ciri.wecon.04)
plot.out.ciri.wopol.02 <- moreprep(plot.out.ciri.wopol.02)
plot.out.ciri.wopol.04 <- moreprep(plot.out.ciri.wopol.04)
plot.out.ciri.wosoc.02 <- moreprep(plot.out.ciri.wosoc.02)
plot.out.ciri.wosoc.04 <- moreprep(plot.out.ciri.wosoc.04)
plot.out.Caprioli_07 <- moreprep(plot.out.Caprioli_07)
plot.out.McDermott_07 <- moreprep(plot.out.McDermott_07)
#out.Clan.11.ols <- moreprep(out.Clan.11.ols)
plot.out.Brinton.10 <- moreprep(plot.out.Brinton.10)
plot.out.McDermott.11 <- moreprep(plot.out.McDermott.11)
plot.out.GII.ols <- moreprepOLS(plot.out.GII.ols)
plot.out.Clan.11 <- moreprep(plot.out.Clan.11)

# pch = 21
# col = black

setwd("../Drafts/Charts")

#pdf(file="ISQFINALMainBodyFigure5.pdf",height=8,width=7)

tiff(file="ISQFINALMainBodyFigure5.tiff",
	res=300,height=8,width=7,units="in")
	
par(mfrow=c(6,2),oma=c(0,0,3,1),mar=c(3,3,1,1))

plot(1,
     pch="",
     ylim=c(1.15*min(plot.out.ciri.wecon.02$Y0),
            abs(1.15*max(plot.out.ciri.wecon.02$Y1))), 
     xlim=c(.5,max(plot.out.ciri.wecon.02$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" CIRI WECON 2002",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.ciri.wecon.02$Y0,
         y1=plot.out.ciri.wecon.02$Y1,
         x0=plot.out.ciri.wecon.02$Index,
         x1=plot.out.ciri.wecon.02$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.ciri.wecon.02$Index,
       y=plot.out.ciri.wecon.02$Estimate,
       col="black",
       pch=21,
      # bg=ifelse(plot.out.ciri.wecon.02$Significance == 0,"white","black"),
      bg=plot.out.ciri.wecon.02$Significance,
       cex=1.5)
axis(side=1,at=plot.out.ciri.wecon.02$Index,
     labels=plot.out.ciri.wecon.02[1:nrow(plot.out.ciri.wecon.02),]$Model,
     cex.axis=1.05, tick=TRUE)


plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.ciri.wecon.04$Y0),
            abs(1.15*max(plot.out.ciri.wecon.04$Y1))), 
     xlim=c(.5,max(plot.out.ciri.wecon.04$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" CIRI WECON 2004",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.ciri.wecon.04$Y0,
         y1=plot.out.ciri.wecon.04$Y1,
         x0=plot.out.ciri.wecon.04$Index,
         x1=plot.out.ciri.wecon.04$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.ciri.wecon.04$Index,
       y=plot.out.ciri.wecon.04$Estimate,
       col="black",
       pch=21,
 #      bg=ifelse(plot.out.ciri.wecon.04$Significance == 0,"white","black"),
        bg=plot.out.ciri.wecon.04$Significance,
 
 cex=1.5)
axis(side=1,at=plot.out.ciri.wecon.04$Index,
     labels=plot.out.ciri.wecon.04[1:nrow(plot.out.ciri.wecon.04),]$Model,
     cex.axis=1.05, tick=TRUE)



plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.ciri.wopol.02$Y0),
            abs(1.15*max(plot.out.ciri.wopol.02$Y1))), 
     xlim=c(.5,max(plot.out.ciri.wopol.02$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" CIRI WOPOL 2002",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.ciri.wopol.02$Y0,
         y1=plot.out.ciri.wopol.02$Y1,
         x0=plot.out.ciri.wopol.02$Index,
         x1=plot.out.ciri.wopol.02$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.ciri.wopol.02$Index,
       y=plot.out.ciri.wopol.02$Estimate,
       col="black",
       pch=21,
       bg=plot.out.ciri.wopol.02$Significance,
       cex=1.5)
axis(side=1,at=plot.out.ciri.wopol.02$Index,
     labels=plot.out.ciri.wopol.02[1:nrow(plot.out.ciri.wopol.02),]$Model,
     cex.axis=1.05, tick=TRUE)


plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.ciri.wopol.04$Y0),
            abs(1.15*max(plot.out.ciri.wopol.04$Y1))), 
     xlim=c(.5,max(plot.out.ciri.wopol.04$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" CIRI WOPOL 2004",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.ciri.wopol.04$Y0,
         y1=plot.out.ciri.wopol.04$Y1,
         x0=plot.out.ciri.wopol.04$Index,
         x1=plot.out.ciri.wopol.04$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.ciri.wopol.04$Index,
       y=plot.out.ciri.wopol.04$Estimate,
       col="black",
       pch=21,
#       bg=ifelse(plot.out.ciri.wopol.04$Significance == 0,"white","black"),
        bg=plot.out.ciri.wopol.04$Significance,
        cex=1.5)
axis(side=1,at=plot.out.ciri.wopol.04$Index,
     labels=plot.out.ciri.wopol.04[1:nrow(plot.out.ciri.wopol.04),]$Model,
     cex.axis=1.05, tick=TRUE)



plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.ciri.wosoc.02$Y0),
            abs(1.15*max(plot.out.ciri.wosoc.02$Y1))), 
     xlim=c(.5,max(plot.out.ciri.wosoc.02$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" CIRI WOSOC 2002",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.ciri.wosoc.02$Y0,
         y1=plot.out.ciri.wosoc.02$Y1,
         x0=plot.out.ciri.wosoc.02$Index,
         x1=plot.out.ciri.wosoc.02$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.ciri.wosoc.02$Index,
       y=plot.out.ciri.wosoc.02$Estimate,
       col="black",
       pch=21,
 #      bg=ifelse(plot.out.ciri.wosoc.02$Significance == 0,"white","black"),
        bg=plot.out.ciri.wosoc.02$Significance,
      cex=1.5)
axis(side=1,at=plot.out.ciri.wosoc.02$Index,
     labels=plot.out.ciri.wosoc.02[1:nrow(plot.out.ciri.wosoc.02),]$Model,
     cex.axis=1.05, tick=TRUE)


plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.ciri.wosoc.04$Y0),
            abs(1.15*max(plot.out.ciri.wosoc.04$Y1))), 
     xlim=c(.5,max(plot.out.ciri.wosoc.04$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" CIRI WOSOC 04",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.ciri.wosoc.04$Y0,
         y1=plot.out.ciri.wosoc.04$Y1,
         x0=plot.out.ciri.wosoc.04$Index,
         x1=plot.out.ciri.wosoc.04$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.ciri.wosoc.04$Index,
       y=plot.out.ciri.wosoc.04$Estimate,
       col="black",
       pch=21,
       #        bg=ifelse(plot.out.ciri.wosoc.04$Significance == 0,"white","black"),
       bg=plot.out.ciri.wosoc.04$Significance,
       
       cex=1.5)
axis(side=1,at=plot.out.ciri.wosoc.04$Index,
     labels=plot.out.ciri.wosoc.04[1:nrow(plot.out.ciri.wosoc.04),]$Model,
     cex.axis=1.05, tick=TRUE)


plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.Caprioli_07$Y0),
            1.15*max(plot.out.Caprioli_07$Y1)), 
     xlim=c(.5,max(plot.out.Caprioli_07$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" Caprioli 2007",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.Caprioli_07$Y0,
         y1=plot.out.Caprioli_07$Y1,
         x0=plot.out.Caprioli_07$Index,
         x1=plot.out.Caprioli_07$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.Caprioli_07$Index,
       y=plot.out.Caprioli_07$Estimate,
       col="black",
       pch=21,
#       bg=ifelse(plot.out.Caprioli_07$Significance == 0,"white","black"),
        bg=plot.out.Caprioli_07$Significance,
        cex=1.5)
axis(side=1,at=plot.out.Caprioli_07$Index,
     labels=plot.out.Caprioli_07[1:nrow(plot.out.Caprioli_07),]$Model,
     cex.axis=1.05, tick=TRUE)


plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.McDermott_07$Y0),
            -1.15*max(plot.out.McDermott_07$Y1)), 
     xlim=c(.5,max(plot.out.McDermott_07$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" McDermott 2007",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.McDermott_07$Y0,
         y1=plot.out.McDermott_07$Y1,
         x0=plot.out.McDermott_07$Index,
         x1=plot.out.McDermott_07$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.McDermott_07$Index,
       y=plot.out.McDermott_07$Estimate,
       col="black",
       pch=21,
#       bg=ifelse(plot.out.McDermott_07$Significance == 0,"white","black"),
        bg=plot.out.McDermott_07$Significance,
        cex=1.5)
axis(side=1,at=plot.out.McDermott_07$Index,
     labels=plot.out.McDermott_07[1:nrow(plot.out.McDermott_07),]$Model,
     cex.axis=1.05, tick=TRUE)



plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.Brinton.10$Y0),
            -1.15*max(plot.out.Brinton.10$Y1)), 
     xlim=c(.5,max(plot.out.Brinton.10$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" Brinton's Discrepancy Measure",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.Brinton.10$Y0,
         y1=plot.out.Brinton.10$Y1,
         x0=plot.out.Brinton.10$Index,
         x1=plot.out.Brinton.10$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.Brinton.10$Index,
       y=plot.out.Brinton.10$Estimate,
       col="black",
       pch=21,
#       bg=ifelse(plot.out.Brinton.10$Significance == 0,"white","black"),
       bg=plot.out.Brinton.10$Significance,
       cex=1.5)
axis(side=1,at=plot.out.Brinton.10$Index,
     labels=plot.out.Brinton.10[1:nrow(plot.out.Brinton.10),]$Model,
     cex.axis=1.05, tick=TRUE)




plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.McDermott.11$Y0),
            -1.15*max(plot.out.McDermott.11$Y1)), 
     xlim=c(.5,max(plot.out.McDermott.11$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main=" McDermott 2011",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.McDermott.11$Y0,
         y1=plot.out.McDermott.11$Y1,
         x0=plot.out.McDermott.11$Index,
         x1=plot.out.McDermott.11$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.McDermott.11$Index,
       y=plot.out.McDermott.11$Estimate,
       col="black",
       pch=21,
#       bg=ifelse(plot.out.McDermott.11$Significance == 0,"white","black"),
        bg=plot.out.McDermott.11$Significance,
        cex=1.5)
axis(side=1,at=plot.out.McDermott.11$Index,
     labels=plot.out.McDermott.11[1:nrow(plot.out.McDermott.11),]$Model,
     cex.axis=1.05, tick=TRUE)



plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.GII.ols$Y0),
           -1.15*max(plot.out.GII.ols$Y1)), 
     xlim=c(.5,max(plot.out.GII.ols$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="UN Gender Inequality Index (OLS Coefficents)",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.GII.ols$Y0,
         y1=plot.out.GII.ols$Y1,
         x0=plot.out.GII.ols$Index,
         x1=plot.out.GII.ols$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.GII.ols$Index,
       y=plot.out.GII.ols$Estimate,
       col="black",
       pch=21,
       #       bg=ifelse(plot.out.GII.ols$Significance == 0,"white","black"),
       bg=plot.out.GII.ols$Significance,
       cex=1.5)
axis(side=1,at=plot.out.GII.ols$Index,
     labels=plot.out.GII.ols[1:nrow(plot.out.GII.ols),]$Model,
     cex.axis=1.05, tick=TRUE)




plot(2,
     pch="",
     ylim=c(1.15*min(plot.out.Clan.11$Y0),
            -1.15*max(plot.out.Clan.11$Y1)), 
     xlim=c(.5,max(plot.out.Clan.11$Index)+.5), 
     xaxt="n",
     xlab="",
     ylab="Coefficient Estimate",
     main="Clan Governance",
     cex.main=1.25,
     adj=1)
abline(h=0,lty="dotdash")

segments(y0=plot.out.Clan.11$Y0,
         y1=plot.out.Clan.11$Y1,
         x0=plot.out.Clan.11$Index,
         x1=plot.out.Clan.11$Index,
         col="black",
         lwd=1,
         lty=1)
points(x=plot.out.Clan.11$Index,
       y=plot.out.Clan.11$Estimate,
       col="black",
       pch=21,
       #       bg=ifelse(plot.out.Clan.11$Significance == 0,"white","black"),
       bg=plot.out.Clan.11$Significance,
       cex=1.5)
axis(side=1,at=plot.out.Clan.11$Index,
     labels=plot.out.Clan.11[1:nrow(plot.out.Clan.11),]$Model,
     cex.axis=1.05, tick=TRUE)



# Title for whole thing
mtext("Oil Income And Indices of Women's Rights", side=3, line=1, outer=TRUE, cex=1.15, font=2)

dev.off()
		


